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j /' j I Transverse magnetic (TM) scattering of an electromagnetic wave from a periodic di- 

electric diffraction grating can mathematically be described by a volume integral equation. 
This volume integral equation, however, in general fails to feature a weakly singular integral 
operator. Nevertheless, after a suitable periodization, the involved integral operator can be 
efficiently evaluated on trigonometric polynomials using the fast Fourier transform (FFT) 
and iterative methods can be used to solve the integral equation. Using Fredholm theory, 

f^ f [ we prove that a trigonometric Galerkin discretization applied to the periodized integral 

(-H ' equation converges with optimal order to the solution of the scattering problem. The main 

advantage of this FFT-based discretization scheme is that the resulting numerical method 
is particularly easy to implement, avoiding for instance the need to evaluate quasiperiodic 
Green's functions. 
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Abstract 



Periodic dielectric structures are important ingredients for modern optical technologies, serving 



> . 1 Introduction 

in 

(N 

O , 

^T) , as beam splitters, lenses, monochromators, and spectrometers. Simulation of electromagnetic 

*y^ ' fields in such periodic structures is a challenging task, since the wave field oscillates in an un- 

^^ ■ bounded domain, since the quasi-periodicity needs to be taken into account, and since evanescent 

^^ i waves arise around the structure. Hence, it might be difficult to use, e.g., a standard finite el- 

ement software for the simulation of wave fields in such structures. For this reason, this paper 
presents a simple-to-implement volume integral equation solver for this simulation task. 
/\^ , We consider scattering of time-harmonic electromagnetic waves from diffraction gratings, 

3 I three dimensional dielectrics that are periodic in one spatial direction and invariant in a second, 

orthogonal, direction (compare Figure [T]). If the incident wave is a transverse-magnetic (TM) 
wave, the electromagnetic field can be described by the scalar equation 

dW{aVu) + k^u = 0, (1) 
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with wave number k > 0, see, e.g., Nedelec ('01)| . The real material parameter a is in this paper 
assmned to be scalar, positive and possibly discontinuous. This periodic scattering problem can 
be equivalently reformulated as a volume integral equation that is formally of the second kind. 
However, since the coefficient a in ([1]) appears in the highest-order term, the integral operator 
of this volume integral equation fails to be compact unless a is not globally smooth (compare, 
for the case of Maxwell's equations, Colton and Kress ('92), Chapter 9.2]). The aim of this 



paper is to analyze the convergence of a trigonometric Galerkin discretization of this volume 
integral equation for discontinuous material parameter a. This analysis will be partly based 
on (purely analytic) results from the paper [Lechleiter and Nguyen ('12)| . Here, we adapt the 
volume integral equations corresponding to ([T]) such that they can be numerically treated via 
an FFT-based approach. This resulting numerical scheme can be rigorously shown to be (quasi- 
optimally) convergent. We provide fully discrete formulas for the implementation of the scheme 
together with computational examples. 




Figure 1: Sketch of the diffraction grating under consideration. 



It might seem inappropriate to consider the TM mode equation ([T]) , since the corresponding 
transverse electric (TE) mode yields the well-known Lippmann-Schwinger integral equation that 
features a weakly singular integral operator. Indeed, the numerical scheme from Vainikko ('00)| 
for the Lippmann-Schwinger equation inspired the scheme we develop here. However, if ma- 
terials feature both dielectric and magnetic contrast then highest-order coefficients cannot be 
avoided even in the TE or TM mode problems. Note that it would not be too difficult to con- 
struct numerical schemes for the simulation of such materials by combining the one from this 
paper with, e.g., schemes developed earlier for the Lippmann-Schwinger equation. 

Volume integral equations are a standard numerical tool in the engineering commu- 
nity to solve scattering problems numerically, see, e.g., Richmond ('65) , Richmond ('66) 



Zwamborn and van den Berg ('92) , Kottmann and Martin ('00) , Ewe et al. ('07) . The linear 



system resulting from the discretization of the integral operator (usually done by collocation or 
finite element methods) is large and dense. Fortunately, the convolution structure of the integral 
operator allows to compute matrix-vector multiplications by the EFT in an order-optimal way 
(up to logarithmic terms), see, e.g., Zwamborn and van den Berg ('92) Rahola ('96)| , at least 
if the discretization respects this convolution structure. This partly explains the success of such 
methods in applications. However, the discretization of the integral operator itself is at least 
in some works done in a mathematically crude way and a rigorous convergence analysis for the 
different discretization techniques is usually missing. 

Despite their relevance in applications, volume integral equations featuring strongly sin- 
gular integral operators (i.e., integral operators that fail to be weakly singular) are a recent 



analytic research subject in mathematics, see, e.g., [Potthast ('99) , Kirsch and Lechleiter ('09) 



Costabel et al. ('10) , Costabel et al. ('12) Lechleiter and Nguyen ('12) . In particular, the nu 



merical analysis of practically feasible discretization methods based on these equations seems 
to be in a somewhat premature stage. Of course, one reason for this phenomenon is that 
for many relevant material configurations, the need for discretizing a strongly singular volume 
integral equation can be avoided. For example, whenever material parameters are piecewise 
constant, boundary integral equations are a powerful alternative to the volumetric approach, 
see, e.g., Otani and Nishimura ('09)| for a recent reference dealing with a periodic scattering 



problem. If the material parameters fail to be piecewise constant, an important approach to 
avoid the discretization of strongly singular integral operators is to combine volume and sur- 
face integral operators. For the full Maxwell's equations in free space, the analytic equivalence 
of both the volume integral equation and the coupled system of weakly singular volume and 
surface integral operators has been worked out in detail in Costabel et al. ('10)|. 



However, whenever using (possibly coupled) boundary integral equations one usually 
needs to be able to rapidly and accurately evaluate the underlying Green's function. It 
is well-known that this is a non-trivial task for (quasi-)periodic Green's functions, see. 



e.g. Linton ('98) , becoming even more difficult if additionally multi-pole expansions are used 
as in Otani and Nishimura ('09)]. The numerical scheme presented here does not require to 



evaluate Green's functions and it is in principle applicable to arbitrary varying material pa- 
rameters. However, the scheme explicitly requires the (two-dimensional) Fourier coefficients of 
the material parameter. According to our experience, the accuracy of computational results 
improves considerably if these coefficients can be computed analytically, or at least be reduced 
to some semi-analytic form that can easily be treated numerically with high accuracy. The 
latter is for instance the case for piecewise polynomial or trigonometric material parameters, as 
we illustrate through examples in the last section. 

Our numerical analysis of a trigonometric Galerkin discretization applied to the 
volume integral equation relies in parts on Garding inequalities that we proved 
in (Lechleiter and Nguyen ('12)] . Of course, these inequalities would in principle directly jus- 
tify any Galerkin discretization of the integral equation. However, such a discretization does 
generally not profit from the above-described advantages arising from the convolution structure 
of the integral operator, the related diagonalization of the operator on trigonometric polynomi- 
als, and the possibility of rapidly evaluating the integral operator using the FFT. Additionally, 
when discretizing the integral operator using finite elements, the strong singularity of the kernel 
makes the computation of the diagonal of the system matrix challenging, see Kone ('10)]. To 



this end, we first periodize the integral operator before discretizing, using a technique that was 
(up to a smoothing procedure) analogously used in [Vainikko ('00)| . The periodized operator 
is then easily evaluated spectrally, since one can (almost) explicitly compute its Fourier coeffi- 
cients (see ([23]) ) . Due to the lack of compactness of the integral operator it seems difficult to 
analyze collocation discretizations as it was originally done in [Vainikko ('00)| . However, it is 
still possible to fully analyze a Galerkin discretization (see Proposition 14. 1 p . 

In essence, the advantage of this trigonometric Galerkin discretization is that it is par- 
ticularly simple to implement - the core of our implementation takes less than 70 lines in 
MATLAB - and that the linear system can be evaluated at FFT speed. By using rela- 



tively simple parallelization techniques on modern multi-core processors this allows to eval- 
uate the integral operator rapidly (MATLAB even automatically uses parallelized FFTW rou- 
tines Frigo and Johnson ('05)| ). Additionally, the FFT-based method requires no evaluation of 
the quasiperiodic Green's function or of its partial derivatives. Due to the slow convergence of 
standard expressions of this Green's function, sophisticated techniques like Ewald summation 
need to be used to accurately evaluate them. Of course, the price to pay for these advantages 
is that the convergence order of this FFT-based method is low if the medium has jumps, due to 
the use of global trigonometric basis functions (otherwise the method is high-order convergent). 
Nevertheless, if one is merely interested in obtaining a moderately accurate solution without 
investing much implementation work, we are convinced that the method presented here is an in- 
teresting simulation technique. This technique could be further improved by using non-uniform 
FFTs that allow some refinement of the underlying grid of the FFT close to edges of the struc- 
ture, for instance. See, e.g., Nie et al. ('05) Zhang and Liu ('02)| for references on non-uniform 



FFTs and their use to solve volume integral equations. 

The rest of this paper is organized as follows: In Section [2] we briefly recall the volumetric 
integral equation for the direct scattering problem and the corresponding Garding inequality 
from [Lechleiter and Nguyen ('12)| . In Section [3] we periodize the volume integral equation such 
that it is suitable for a fast FFT-based discretization on biperiodic trigonometric polynomials. 
We also prove the necessary Garding inequalities for the periodized system (see Theorem 13. 5p . 
These inequalities are naturally the basis for quasi-optimal error estimates for the trigonometric 
Galerkin discretization in Section HI Finally, Section O contains several illustrative numerical 
examples. 

Notation: L^-based Sobolev spaces on a domain D are denoted as H^{D), s G M, and 
C^'^{D) is the usual space of Lipschitz continuous functions that possess Lipschitz continuous 
partial derivatives up to order m. Further, Hf^^(D) = {v £ H^(B) for all open balls B C D}. 
The trace of a function u on the boundary dD from the outside and from the inside of D is 
denoted as 7cxt(^) and 7int(^i), respectively. The jump of u across dD is [u]q£) = 7cxt(^)— 7int('u). 
If the exterior and the interior trace of a function u coincide, then we simply write 7(m) for the 
trace. 

2 Problem Setting and Known Results 

Propagation of time-harmonic electromagnetic waves in an inhomogeneous, isotropic, and loss- 
less medium is described by the Maxwell's equations for the electric and magnetic fields E and 
H, respectively, cnilH + iueE = and curl£^ — iuj^qH = in M^. Here, oj > denotes the 
frequency, e is the positive electric permittivity and /Uq is the (constant and positive) magnetic 
permeability. We assume in this paper that the scalar function e is independent of the third 
variable X3, and 27r-periodic in the first variable xi. Further, we suppose that e equals a constant 
eo > outside the grating structure. 

If an incident electromagnetic plane wave independent of the third variable X3 illuminates 
the grating, then the Maxwell's equations for the total wave field decouple into two scalar partial 
differential equations. In particular, the third component ^^3 of the magnetic field satisfies 

div iyE^^Vu) + k^u = with e^ := e/eq and k := w-y/eo/Uo > 0, (2) 



together with jump conditions on interfaces where the refractive index e~^ jumps: u and 
e~^du/dv are continuous across such interfaces. Note that Er is 27r-periodic in xi. We as- 
sume that the contrast q := e^^ — 1 has support in {\x2\ < p} for some constant p > 0. 

Consider now a plane incident wave u^{x) = exp{ikx ■ d) = exp{ik{xidi + X2d2)) where 
|(i| = 1 and ^2 7^ 0. When n* illuminates the diffraction grating there arises a scattered field u^ 
such that the total field u = u^ -\- u^ satisfies ([2]), that is, the scattered field satisfies 

div {e-^Vu") + k^u' = -div (gVn*) in M^. (3) 

Note that u* is a-quasi-periodic with respect to xi, 

u\xi + 27r,X2) = e^'^^'^u\xi,X2) for a := kdi. 

Since Er is periodic we seek for a scattered field that is a-quasi-periodic in xi, too. For unique- 
ness of solution we require that u^ above (below) the dielectric structure can be represented 
by a uniformly converging Rayleigh series consisting of upwards (downwards) propagating or 
evanescent plane waves, 



u^(^x) = ^ufe"'^^'^"'^^^^^f\ X2^±p, aj:=j + a, /3,- := ^fc^ - a|. (4) 



The square root used to define 



Pi = V ^ - "^ '■= U^H I2M/2 ,2 / J ' J e 



_i(a2-fc2)l/2^ J,2^^.^ 



is chosen such that Im(/3j) > always. Further, the so-called Rayleigh coefficients u- of the 
scattered wave in Q have explicit representations. 

Note that we call a solution to the Helmholtz equation radiating if it satisfies dU . 

By Ga we denote the Green's function to the a-quasi-periodic Helmholtz equation in M^, 



see, e.g., Linton ('98), Eq. (2.13)]. In this paper, 

we always suppose that k / a^ for all j G Z, (5) 



which implies that this Green's function has the series representation 

Gq,(x) := — ^ — exp(iajXi -hi^j|x2|) for x = f M , x / f ^ j,mG 



(6) 



Note that ([5]) implies that all the /3j = (A:^ — a^)^" are non-zero, and that the Green's function 
is well-defined, see again Linton ('98)]. 



Remark 2.1 (Failure at Wood's anomalies). The phenomenon that condition ^ fails to hold 
for some k > is called a Wood's anomaly, see, e.g., Barnett and Greengard ('11)^. At 



a 



Wood's anomaly, the representation ([6]) is obviously not well-defined. Image-like representa- 



tions of the Green's function would (at least formally) he well-defined, see, e.g., Linton ('98) 
Eq. (2. 7)] for an example. Hence, it might seem as if there was a chance that the method pre- 
sented in this paper works at Wood's anomalies. However, by carefully checking Lemma I3.il 
below one notes that this is not the case since certain Fourier coefficients (denoted by ICp later 
on) are not well-defined at Wood's anomalies. 

We introduce the strip Q := (— 7r,7r) x M and set 

nn.= {-iT,TT)x {-R,R) forR>0. 

Moreover, we set H^{Q,ji) := {u G H^{Q,ji) : u = U\aj^ for some a-quasi-periodic U S 
-?/jqj,(M^)} for ^ E N and i? > 0, and H^{Q) is defined analogously. For any Lipschitz do- 
main D (see [McLean ('00)| for a definition), the space L^(L',C^) contains all square integrable 
functions with values in C'^ (complex column vectors with two components). 

Lemma 2.2 (Lemmas 5 and 6 in [Lechleiter and Nguyen ('12)| ). If D C ^ is a Lipschitz 
domain, then the volume potential 



(vm^) 



/ Ga{x-y)f{y)dy , x £ Qr, 
Jd 



is bounded from L^{D) into H^(Q,ji) for all R > 0. For g G L^{D, C^) the potential w = div Vg 
belongs to H\{Q.ji) for all R > 0. It is the unique radiating weak solution to Aw-'rk'^w = — div^f 
in Q, that is, it satisfies the Rayleigh expansion condition Q, and 

I (V. . V^ - e.-v) d. = - / . . VWd. for aU v e H},ia) «ntk cor.,act >u^n. (7) 

Let us now come back to the differential equation Q for the scattered field u^ . Recall that 
we assumed that the contrast q = £~^ — 1 has support in {\x2\ < p} for some p > 0. We denote 
this support (restricted to one period — vr < xi < vr) by 

D = supp(g) 

and suppose from now on that D is a Lipschitz domain. If we choose R > p, then D C ^r. 
Moreover, by setting / = gVn* in ^ the variational formulation of ^ reads 

/ (Vn' • VU - k'^u'v) dx =- [ {qVu' + /) • Vvdx (8) 

Jn Jd 

for all V G H^[Q) with compact support in Q,. From Lemma 12.21 we know that the radiating 
solution to this problem is given by u^ = div V{qVu^ + /)• If ^^ define the bounded linear 
operator 

L: L\D,C^)^Hi{D), f^diwVf, 

6 



then the scattered field n^, solution to ([3D, hence solves the volume integral equation 

u^ - L{qVu^) = L{f) mHl[D) (9) 

for / = qVu^. The operator on the left of the last equation satisfies a Garding inequality. 



Theorem 2.3 (Theorem 16 in Lechleiter and Nguyen ('12)|). Assume that q > qq > in D, 



that y^ G C^'^(Z)), and that D is of class C^'^. There exists a compact operator K on H^{D) 
such that 

Re{v - L{qVv),v)Hi(^D) > 11^11^1 (d) -^e{Kv, f)Hi{D), ve H^{D). 

For a real- valued contrast, uniqueness of solution of the scattering problem dSHH), or equiv- 
alently of the integral equation ([U]), does in general only hold for all but a discrete set 
of positive wave numbers. Uniqueness results either require (partially) absorbing materi- 
als or non-trapping coniditions on the material; examples of such conditions are given, e.g., 

Bonnet-Ben Dhia and Starling ('94) .Elschner and Schmidt ('98)|. 



m 



Remark 2.4 (Assumption on uniqueness of solution). In the rest of the paper, we always 
suppose that uniqueness of solution to (Oj^P holds. 

We restrict our theoretical analysis to real and positive contrasts, since Garding 
inequalities corresponding to complex-valued or negative contrasts are more involved, 
see Lechleiter and Nguyen ('12)| . Treating these cases would increase technicalities without 
adding new ideas to the text. 

3 Periodization of the Integral Equation 

In this section we periodize the volume integral equation Q and show the equivalence of the 
periodized equation and the original one. The purpose of this periodization is that the resulting 
integral operator is, roughly speaking, diagonalized by trigonometric polynomials. This allows to 
use fast FFT-based schemes to discretize the periodized operator and iterative schemes to solve 
the discrete system. We also prove Garding inequalities for the periodized integral equation, 
which turns out to be involved. However, these estimates are crucial to establish convergence 
of the discrete schemes later on. 

Let us again emphasize that we assume in all the paper that the non-resonance condition ([5]) 
is satisfied, which excludes Wood's anomalies. 

Since we are interested in spectral schemes we define a periodized Green's function, firstly 
setting 

ICp{x) := Gaix), x = (2;i,X2)'^ G M X (-p,p), X / (27rm,0)^ for m G Z, (10) 

and secondly extending K,p{x) 2p-periodically in X2 to R^. The trigonometric polynomials 



fj{x)-=—r^e^V\^{h + a)xi+\ X2), j = {ji,J2) GZ, (11) 



are orthonormal in L'^{Q.p). They differ from the usual Fourier basis only by a phase factor 
exp(iaxi), and hence also form a basis of L'^{VLp). For / S L'^{Vtp) and j = {ji,J2)^ £ Z^, 
f[j) := J"^ fTpjdx are the Fourier coefficients of /. For < s < oo we define a fractional 
Sobolev space -ffpgj.(ilp) as the subspace of functions in L^(r2p) such that 



A^p) •" 



y:(i+ 



v(i)r<oo. 



jez2 



/orfcV(ii+«)'-(^) 
eZse, 



J 



G 



It is well-known that for integer values of s, these spaces correspond to spaces of a-quasi-periodic 
functions that are s times weakly differentiable, and that the above norm is then equivalent to 
the usual integral norms. 

Lemma 3.1 (Theorem 2 in Lechleiter and Nguyen ('12)| ). The Fourier coefficients of the ker- 
nel ICp from (fTUj) are given by 

! l cos{J2Tt) exp(i/3j^ p)-l ^ ,9 / /• , x9 / ioTT\2 
J_ (p^3/2 
4J2 Vtt/ 

The convolution operator Kp, defined by {Kpf)[x) = J^ K.p{x—y)f{y) dy for x G ftp, is bounded 
from L'^{Q.p) into Hp^j.{flp). 

The periodized kernel ICp from pO|) is not smooth at the boundaries {x2 = ip}- To prove 
Garding inequalities for the periodized integral equation, we additionally need to smoothen this 
kernel at {x2 = ip} and, to this end, introduce a suitable cut-off function. For R > 2p we 
choose a 2i?-periodic function x £ C^(]R) that satisfies < x < 1 and x(^2) = 1 for |x2| < 2p. 
Moreover, we assume that x(i?) vanishes up to order three, x i^) = for j = 1, 2, 3 (compare 
Figure [2D. 
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Figure 2: (a) The 2i?-periodic function x equals to one for \x2\ < 2p, and it vanishes at ±i? up 
to order three. In this sketch, p = \ and i? = 4. (b) The support of the contrast (shaded) is 
included in i7p = {|x2| < p}, and R > 2p. 



Let us define a smoothed kernel /Csm by 



ICsm{x) = x{x2)ICr{x) for x G M^, x ^ I ^^; j ' ^ ^ 



/^27rii\ 



(12) 



where ICr is the kernel from (jlOp . Note that ICsm is a-quasi-periodic in xi, 2i?-periodic in X2, 
and a smooth function on its domain of definition (that is, away from the singularity). 



Lemma 3.2. The integral operator Lpcr : L^{i^R,C^) — )■ H1^^^{Qr) defined by 

Lperf ■■= div / ICsm{- " y)/(y) dy 

is bounded. 

Proof. We split the integral operator in two parts, 

Lpcrf = div / /Csm(- - y)fiy) dy = div / x(- - y2)K,R{- - y)f{y) dy 



= div/ ICR{--y)f{y)dy +div [xi- - y2) - MICr{- - y)f{y) dy . 

By Theorem 13. H the integral operator with the kernel ICr is bounded from L^(J7/j,C^) into 
H^{Ur). Further, the definition of x shows that x(x2 — 2/2) = 1 for \x2 — 2/2! < 2/9. The kernel 
(x — 1)/C_R is hence smooth in Qr, and the corresponding integral operator is compact from 
L^(0/{,C^) into H^(^[}r). Hence, Lpcr is bounded from L'^{Qr,C'^) into H^{Ur). Periodicity of 
the kernel /Cgm in the second component of its argument finally implies that Lper/ belongs to 

hI,,{Qr) c hHQr). a 

Let us now consider the periodized integral equation 

u - Lp,r{qVu) = Lp,,{f) in i7^er(^ij), (13) 

where, for simplicity, we call the unknown function u. 

Theorem 3.3. (a) If f € L'^{D,<C?), then Lpcr(/) equals L{f) in Qp. 

(b) Equation ([9]) is uniquely solvable in H^{D) for any right-hand side f E L^(Z),C^) if and 
only if (|13p is uniquely solvable in H^^^{Qr) for any right-hand side f S L'^{D,'C?). 

(c) If q £ C^'^{D) and if f = qVu^ for a smooth a-quasi-periodic function u*, then any 
solution to (fT3]l belongs to ifpgr(^fl) for any s < 3/2. 

Proof, (a) For all x and y € Qr_ such that \x2 — y2\ < 2p it holds that /Csm(a; — y) = x(^2 — 
y2)K,R{x — y) = Ga{x — y). In particular, for x G Jlp and y G D C fip it holds that |x2 — y2| < 2p. 
Consequently, 

(-^pcr(/))(x) = div / /Csm(2; - y)f{y) dy 
Jur 



= div / Ga{x - y)f{y) dy = {L{f)){x), x G Qp. 
Jd 

(b) Assume that u^ G Il\{D) solves ([9|) for a right-hand side / G I?'{D.,€?) and define 
u G H^^^{Qr) by u = LpcriqVu^ + /). Since u^ solves ([9]), and due to part (a), we find that 
u\d = u^- Hence -Lpcr(gVu) = Lpej-{qVu^) in Hp^j.{Qr), which yields that 

u = Lper((?Vn + /) in F^.^C^ij). (14) 



Now, if / G L^(L',C^) vanishes, then uniqueness of a solution to ^ imphes that u^ G H^(D) 
vanishes, too. Obviously, u = Lper{qVu^) vanishes, and hence (fH|) is uniquely solvable. The 
converse follows directly from (a). 

(c) Assume that u G Hj^^,j.{^ji) solves ()13p for / = qVu^. Part (a) implies that the restriction 
of u to Op solves u — L{qVu) = L^qVu^) in H^^iytp). Hence, Lemma 12.21 implies that u is a 
weak a-quasi-periodic solution to div((l + q)Vu) + k'^u = —di\ {q\Iu^) in Q.p. Transmission 
regularity results imply that u belongs to H^(D) n H'^iVlp \ D), and it is well-known that this 



implies that u G H'^{^p) for s < 3/2 (see, e.g., [Grisvard ('92) Section 1.2]). The function u is 



even smooth in Qji \ flp-g: Recall that p was chosen such that D C flp. Hence, there is e > 
such that D C Qp-2e, and 



u(x) = Lpcr(gV(u + u*))(x) = div / ICsmix - y)q{y)V {u{y) + u' (y)) dy , x e^B,\np^s 

Jd 

shows that the restriction of u to VLr \ ^p-e is a smooth a-quasi-periodic function, since the 
kernel of the above integral operator is smooth. D 

Next we prove that the operator / — Lpci-{qV-) from (jlSp satisfies a Garding inequality in 
Hp^j.(i^fi). First, we announce a simple lemma that is useful in the next proof. 

Lemma 3.4. Suppose that X and Y are Hilbert spaces. Let A, B he hounded linear operators 
from X into Y and consider the sesquilinear form (n, v) i— )• {Au, Bv)y on X x X . If either A 
or B is compact, then the linear operator Q : X ^ X defined hy {Qu,v)x = {Au,Bv)y for 
u,v £ X is compact, too. 

Theorem 3.5. Assume that yjq G C^'^(D), that q > qo > 0, ctnd that D is of class C'^'^. Then 
there exists C > and a compact operator K on Hp^j.{Q,ji) such that 

Re{v-Lp^r{q'^v),v)Hi^^^n^) > WvWhi^^^q^) -Re{Kv, i;)//i^^(n^), v G F^^,(Jlfl). (15) 

Remark 3.6. The idea of the proof is to split the integrals defining the inner product on the 
left of (|15p into the three integrals on D, Qp\ D, and on Q^ \ Qp. For the term on D one 
exploits the Garding inequalities from Theorem \2.3[ . The terms on rtp\D and on Q,r \ Vtp can 
he shown to he compact or positive perturhations. 

Proof. Let v G Hp^^{Q.ii). First, we split up the integrals arising from the inner product on 
the left of ([15]) into integrals on D, on Vtp \ D, and on rj/j \ Qp. Second, we use the Garding 
inequality from Theorem 12.31 to find that 

Re{v-Lp,,{qVv),v)Hi^^in^) > ll^llii(D) + (^^. ^)hi(Z)) + lbll^i(c«\D) 

- Re [(Lper(gVf ), v)^,^^^^^^^^ + {Lp,,{qVv), f)j:^i (j^^^^)] (16) 

with a compact operator K on H^{D). Further, the evaluation of Lpcr(Q'V-) on il/j \ Op defines 
a compact integral operator mapping H^{D) to H^{flji \ Qp), because the (periodic) kernel of 
this integral operator is smooth. (This argument requires the smooth kernel /Csm introduced in 
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the beginning of this section.) Lemma 13.41 then allows to reformulate the corresponding term 
in ()16p in the way stated in the claim. Unfortunately, the last term in ()16p does not yield a 
compact sesquilinear form and needs a more detailed investigation. 

For X £ il.p\ D and y £ D the kernel lCsm{x — y) equals Gct{x — y), which is a smooth 
function of x G i7p \ D and y £ D. Moreover, AGa{x — y) + k'^Ga{x — y) = for x 7^ y. Since 
VxGaix — y) = —VyGa{x — y), an integration by parts in Qp\ D shows that 

L{qV v) (x) =div / Ga{x - y)q{y)Vv{y)dy 
Jd 

= - VyGa{x-y)-V{qv){y)dy + / VyG^ix - y) ■Vq{y)v{y)dy 
Jd Jd 

= -k^ Ga{x-y)q{y)v{y)dy - L{vVq){x) 
Jd 

fit I T" II 1 

'^-l—\ l-^rA{q){y)l{v){y) ds for x£VLp\D, 

dD ou{y) 

where v is the exterior normal vector to D. The integral operator appearing in the last term of 
the last equation is the double layer potential DL, 

DLW=/ ^%^T^V(y)ds mn\dD. 
JdD ou{y) 

It is well-known that DL defines a bounded operator from HJ (dD) into H^{Q,ji \ D) and into 
H\{D) (see, e.g., [Arens ('10)| ). This implies that the jump of the double- layer potential 

TV' := [DLV'laz) = 7cxt(DL^) - 7int(DLV) 

1/2 
from the outside of D to the inside of D is a bounded operator on Hd {dD). It is well-known 

that in our case T is even a compact operator on HJ {dD), since D is of class C^'^. Additionally, 
the equality 7int(DL'(/') = --0/2 + Tip holds for V G H^'^{dD). 
For V G HI^^I^r), 

- (VL(gVt-), V7;)^2(^^\;d) = {k'VV{qv) + VL{vVq) + VIih{y,,,,{qv)) , Vv) ^.^^^^^-^y (17) 

The mapping properties of V shown in Lemma 12.21 and the smoothness of q imply that v 1— )• 
k'^\/V{qv) + S/L{vVq) is compact from H^^,^{Qji) into L'^{D). To finish the proof we show that 
the last term in (J17p can be written as a sum of a positive and compact term. For simplicity, we 
define w = DL(7int(gw)) and note that —v/2 = [7int(^'^) — ^(7int(Q'w))]/7int(9) on dD. Since it 
plays no role whether the normal derivative dw/du is taken from the inside or from the outside 
of D, we skip writing down the trace operators for the normal derivative. Then 



(VDL((?^), V^)^2(f,^\I5) = / _ Vt/; • VUdx 



np\D 

, 2 /" - T f dw_ f dw _ f dw _ 

k / wvdx — / —-vds + / 7- — vds — / — — uds (18) 

J^f.\D JdD dv 7p^ dx2 7r_p dx2 
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and the above jump relation shows that 

1 f dw_ , f 
— - / ——vas = / 

2 JsD ov Jq 



^■^^j„ _ f 9wjint{w) _ f dwT{jint{qv)) , 



d^ 7int(g) JdD du 7int(g) 



D V// JD (1 JdD du Trntiq) 

, ^^dx+ f {Vq^'.Vw-k^^)wdx- [ gf! ^(WH) d.. 
Id Q Jd Q JdD dv 7int(g) 

Combmmg the last computation with (jlSp shows that 



(VDL(gt>|az)), ^v\^^\D)L'^{n,\D) = ^ / ^— ^^ + ^ / _"^^dx (19) 




JZ) V 9/ J9D 5^^ 7int(g) 

Using Lemma 13. 4| all the terms in the second line of the last equation can be rewrit- 
ten as {Kiv, v)hi m^\ where Ki is a compact operator on H^^j.{^ji). The mapping v i— )• 
fjj\Vw\^/qdx is obviously positive ii q > 0. In consequence, (fT6]l and (fT7|) show that (fT5]l 
holds. D 

4 Discretization of the Periodic Integral Equation 

In this section we firstly consider the discretization of the periodized integral equation (|13p in 
spaces of trigonometric polynomials. If the periodization satisfies certain smoothness conditions 
and if uniqueness of solution holds, convergence theory for the discretization is a consequence 
of the Garding inequalities shown in Theorem 13.51 Secondly we present fully discrete formulas 
for implementing a Galerkin discretization of the Lippmann-Schwinger integral equation (J13p . 



For iV e N we define Z% = {j £ T? : -N/2 < ji,2 < N/2} and Tn = spanjc^j : j £ Z%}, 
where ipj G L?'{VLji) are the a-quasi-periodic basis functions from (jlip . Note that the union 
UatsnT/v is dense in H^^^{Q.r). The orthogonal projection onto T/v is 

Pn : hI,,{VLr) -^ Tn, Pn{v) = Y, HJ)y^j, 

where v{j) denotes as above the jth Fourier coefficient. The next proposition recalls the stan- 
dard convergence result for Galerkin discretizations of equations that satisfy a Garding inequal- 



ity, see, e.g. [Sauter and Schwab ('07) Theorem 4.2.9], combined with the regularity result from 



Theorem IS.Sfc 



Proposition 4.1. Assume that q satisfies the assumptions of Theorem, \3.5\ and that ([9|) is 
uniquely solvable. Then ()13p has a unique solution u € Hp^,j.{flfi), and then there is Nq G N 
such that the finite- dimensional problem to find u^ £ T/v such that 

{uN - Lpcr{q^UN),WN)H^^,{^n) = if, w^Jv)Hi„,{Qfl) for all wn&Tn (20) 
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possesses a unique solution for all N > Nq and f € H^^j.(flfi). In this case 

\\un - u\\jji (n„) < C inf \\wn - u\\ui (nt,) < CN~^\\u\\rri+s,r^ n, < s < 1/2, 

with a constant C independent of N > Nq. 

Remark 4.2. The convergence rate increases to s+l — t if one measures the error in the weaker 
Sobolev norms of Hpgj.{Q,ji), 1/2 < t < 1. This can be shown using adjoint estimates (see, 
e.g. Sauter and Schwab ('07). Section 4-2] for the general technique). However, the (linear) 
rate saturates at t = 1/2, since the integral operator is not bounded on H^^^{Q.fi) for t < 1/2, 
that is, the L'^ -error decays with a linear rate. 

Applying P/v to the infinite-dimensional problem ()13p . and exploiting that Pn commutes 
with the periodic convolution operator Lper, we obtain the discrete problem to find njv G T/v 
such that 

UN - Lpcr(-PAr(q'VuAr)) = Lpcr(-P7v/). (21) 



Fast methods to evaluate the discretized operator in ()2ip exploit that the application of Lpcr 
to a trigonometric polynomial in T/v can be explicitly computed using an a-quasi-periodic 
discrete Fourier transform that we call J-'jy. This transform maps point values of a trigonometric 
polynomial (pj (see (jlip ) to the Fourier coefficients of the polynomial. If /i := {2it/N, AttR/N)'^ 
(a column vector), then 

%(j) = -^^ ^vn{1- h) exp ( - 27ri {h + a, ja)^ • l/N) , j G Z^. 

(Vectors j, / € Z^ are interpreted as a column vectors.) This defines the transform J-jv mapping 
{vn{J ■ h))j(zzl *° (^NU))jez%- The inverse J"^^ is explicitly given by 

vnU ■ h) = -== Y^ vnH) exp (27ri (^i + q, ^2)^ • j/N^ , j G 
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l&% 



Both J^Tv and its inverse are linear operators on C^ = {{cn)n<^i? '■ c„ G C}. The restriction 
operator Rn,m from C^ to C|^, N > M, is defined by Rn,^!^^) = b where b{j) = a{j) 
for j G Z|j. The related extension operator Em,n from C|^ to C^, M < A^, is defined by 
Em,n{o) = b where b{j) = a{j) for j G Z|^ and b{j) = else. 

For the next lemma, we introduce the notation A»B = iAijBij)f'-'^ for the componentwise 
product of two matrices A,B^ (^MxM _ 

Lemma 4.3. The Fourier coefficients of qdiUN, i = 1,2, are given by 

{qdeUN{j))j^Z2^ = R3N,NJ'3N[j'3Ni^'iN,3N{q2N{j))j(.z%) • ^3Ni^N,3N{Wi{j)uN{j))j(zZ%)] 

where wi{j) = i{ji + a) and W2{j) = iJ2'/r/i? for j G Z^. 
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Proof. For un G Tn, j G Z^, and i = 1,2, 

Att R qdiU j^ (j) = 4:71 R / qdiU^'ipj dx = AttR N^ diU]\f{m) / qTpjipmdx (22) 

= (47ri?)i/2 ^ §^^(m)g(j-m). 



-(J2-m2)a;27r/_R] ^^ 



If j € Z^, then the coefficient qdeuj\f{j) merely depends on q{m) for m £ Z^Af- Hence, 
qdiUNiJ) = q2NdeUN{j) for j G Z^. Obviously, q2NdiUN belongs to Tsn- Hence, the Fourier 
coefficients of q2NdeUN are given by T-sn applied to the grid values of this function at j • h, 
j e Zg^. The grid values of B^un are given by T.^p^{EN,3N{diU^{j)j^^2 ), and the grid values 
of q2N can be computed analogously. Finally, taking a partial derivative with respect to xi 
or X2 of u yields a multiplication of the jth Fourier coefficient u{j) by i(ji + a) and ij2'/r/i?, 
respectively. D 

In Lemma 13.11 we computed the Fourier coefficients of the kernel ICji. The kernel /Cgm 
used to define the periodized potential Lpcr is the product of ICr with the smooth function x 
(see (fT2]) ). Hence, the Fourier coefficients of /Csm are convolutions of the JCr{j) with x(j2) = 
(47ri?)-V2 f^^exp{-iJ27rx2/R)x{x2) da^a , 

^sm(j) = /4^oNi/2 X] ^Rih,m2)x{J2 - m2), j G Z^. 

The latter formula can be seen by a computation similar to ()22p . Note that x is a smooth 
function, which means that the Fourier coefficients x in the last formula are rapidly decreasing, 
that is, the truncation the last series converges rapidly to the exact value. The convolution 
structure of Lpcr finally shows that 

(VJ)(i) = (47ri?)i/2/Csm(j)[i(ji+a)/i(i) + ^/2(j)], f=[Q eL\^R,C'). (23) 

The finite-dimensional operator un i— )■ Lpcx{Pn{q^un)) can now be evaluated in 
0(A^log(A^)) operations by combining the formula of Lemma 14.31 with ()23p . The linear sys- 
tem (|21|) can then be solved using iterative methods. Whenever one uses iterative techniques, 
one would of course like to precondition the linear system. The usual multi-grid preconditioning 
technique for integral equations of the second kind (see, e.g., [Vainikko ('00)| ) does not apply 
here, since the integral operator is not compact. For the numerical experiments presented in the 
following section, we simply used the (unpreconditioned) GMRES algorithm from Kelley ('95)|. 
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5 Numerical Examples 



In this section we present several numerical examples that first check the quasi-optimal conver- 
gence rate of the trigonometric Galerkin scheme under investigation from Proposition 14.11 The 
second aim of the examples is to show that the scheme is able to cope with spatially varying 
material configurations and hence can be applied when, e.g., boundary integral techniques are 
blocked. Third, we try to illustrate memory needs and computation times to give an impression 
about the performance of the scheme. Before going into details, let us emphasize again that one 
of the main advantages of the trigonometric Galerkin scheme is its straightforward and rather 
easy implementation, at least if a performant FFT routine and an efficient linear solver are 
already at hand. 

We first present a computational example for a very simple strip-like structure for which 
one can compute the scattered field for an incident plane wave analytically. Checking that 
the numerically computed solution to the scattering problem converges well to the analytically 
known expression will be the first test for the correctness of the code (and a - due to its simplicity 
limited ~ test for the convergence rate, too). Recall that we aim to compute the scattered field 
for an incident field u^{xi,X2) = exp{ik{cos{6)xi — sm{6)x2)) with incident angle 6. For this 
first test we choose k = 7r/2 and 9 = 7r/4 and approximate the solution in T/v where N = 2"^ 
for n = 6, ...,11. For this example, D = (— 7r,7r) x (—0.75,0.75) is a strip, we choose Qpt = 
(— TT, vr) X (—2,2), and the contrast q equals two in D (compare Figure [3l^a)). For this setting 
one can explicitly compute the scattered field and use the analytic expression for comparison. 
In the Figure O^b) we show the relative error between the numerical and the analytical solution 
in the norms H^^^{Qji) where s = 0,0.5, 1. The relative error measured in the norm Hpgj.{Qji) 
fits quite well to the theoretical statement in Proposition 14.11 Furthermore, if one measures 
the relative error in the norm Hp^j.(^[}jj) for s = and s = 0.5 the experiment confirms the 
statement of Remark 14.21 All computations in this and in all following numerical experiments 
were done on a machine with an Intel Xeon 3200 quad core processor and 12 GB memory. 
The trigonometric Galerkin discretization of the volume integral equation was implemented 
in MATLAB, relying on FFTW routines [Frigo and Johnson ('05)] that MATLAB is able to 
execute in parallel. The linear system was solved by the GMRES iteration from Kelley ('95)| . 
The iteration was stopped when the relative residual reduction factor was less than 10~^ (this 
parameter was chosen for all later experiments). Figure [1] shows computation times and the 
number of GMRES iterations for this numerical test. Obviously, the computation time of the 
scheme gets large when N becomes very large due to memory needs. The number of GMRES 
iterations slowly decreases in N from 7 to 5. 



N 



64 128 256 512 1024 2048 



Computation time(s) for strip-structure 
jj GMRES iterations for strip-structure 



0.3 1.1 3.7 21.2 131.7 463.7 
7 6 6 6 6 5 



Table 1: Computation times and number of GMRES iterations for the computation of the 
errors for the simple strip structure shown in Figure [3j The parameter N is the discretization 
parameter of the trigonometric Galerkin scheme. 



15 




i5 

03 

DC 



10" 



10" 



10" 




'□^ 
^■v. 



64 128 256 512 1024 2048 
N 



(a) 



(b) 



Figure 3: (a) Two periods (in the horizontal variable) of the strip structure with contrast equal 
to two. (b) Relative error between the numerically approximated solution and the analytically 
computed reference solution measured in i?pgj,-norm for scattering from the strip shown in (a). 
Circles, kites, triangles correspond to s = 1, s = 0.5 and s = 0, respectively. The continuous line 
and the dotted lines indicate the convergence order 0.5 and 1, respectively. The discretization 
parameter is A^ = 2" for n = 6, ..., 11. 

Of course, the numerical results for the strip structure from the last example merely provide 
a first test that the algorithm computes correct solutions. For further tests and illustrations 
of the algorithm, we consider more complicated structures where the contrast varies smoothly 
within subdomains and jumps across subdomain borders. As we mentioned in the introduction 
and confirmed in Lemma 14.31 it is essential for the Galerkin scheme to have explicit values 
of the Fourier coefficients q2N of the contrast q at hand. In principle, these values could be 
approximated using FFTs. However, we found that whenever one is able to compute these 
Fourier coefficients analytically, this results in considerably more accurate computations. In the 
examples below, we explain case-by-case how to compute these Fourier coefficients for a wide 
class of polynomially or exponentially varying materials. For complicated material shapes, it is 
usually impossible to compute the Fourier coefficients explicitly. Using partial integrations, one 
is however able to come up with semi-analytic expressions that merely require a one-dimensional 
integration of a periodic and piecewise analytic function for evaluation. 

Figure|l]shows the four contrasts 51,2,3,4 that we consider in the experiments below. We start 
now by giving precise definitions of these four contrasts and we compute their Fourier coefficients 
(semi-)explicitly. Afterwards, we present numerical examples for the different structures. We 
would like to point out in advance that for all four examples the domain Jljj will be chosen as 
(— vr, vr) X (—2,2), i.e., R = 2 always. 

The contrast qi plotted in Figure IHa) consists of 27r-periodically aligned kite-shaped 
inclusions with constant material parameter (the contrast equals to two inside the inclu- 
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Figure 4: The four plots show two periods in the horizontal variable of the contrasts ?i,2,3,4 
considered in the numerical experiments below, (a) The piecewise constant kite-shaped contrast 
qi. (b) The piecewise constant contrast (72 is supported in a strip, (c) The contrast (73 varies 
smoothly within a sinusoidally-shaped strip, (d) The contrast q^ varies smoothly within a 
rectangle. 

sion). The boundary of the central inclusion D C (— vr, vr) x (—2,2) is parametrized by 
t I— )■ (1.5cos(t) + cos(2i) — 0.65, sin(i))~'', t G [0, 27r]. The Fourier coefficients {qi{j)}j£z'^ can 
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simply be computed using Green's formula, 



= 2 [ e-'^i^i-'^^^da; =— [ i/2(2;)e-'^i^i-^¥^2ds 
Jd J2vr JdD 

= — I "e-'J'i"iW-'¥'^"^W(1.5sin(t) + 2sin(2t))dt, j2 / 0. 
j2vr7o 

This integral can now be accurately evaluated numerically (we use the fourth-order convergent 
composite Simpson's rule). 

Similar techniques yield the Fourier coefficients of the contrast 

f 1 in Di := (-7r/2,7r/2) x (0,0.75), 

02= { 

[2 in (-vr, vr) x (-0.75, 0.75) \ L»i, 
that is plotted in Figure IDj^b). The Fourier coefficients of (72 can be computed explicitly, 

Jdi Jd 

Both integrals can of course be computed analytically, the first one equals for instance 
4i/(vrjij2) sin(ji7r/2)[l - exp(-3/87riji)] for ji,2 ^ 0. 

The contrast q^, shown in Figuredl^c) is defined as a smooth function on a sinusoidally shaped 
strip D. In detail, 

D = { (xl) £ M."^ ■■ -vr < xi < vr, sin(22;i) < 2x2 - 1 < sin(2xi)} and 

, ^ /e-V3 forx=(^i)GD, 
'^^^"^ = \0 else. 

In this case the Fourier coefficients of the contrast q can be computed semi-analytically using 
Green's formula 



Stt qsij) = [ q{x)e~'^^'''~'^^''''^/'^ dx =- [ e^'^^''^^'-^+'^^^/^^''^ dx 

V3 f ^^^^^^-ij,xi-{l+ij27r/2)x2 ^g 



1 + ij2vr/2 Jq^ 

~^/^ / "" g-ijit-(l+iJ2vr/2)(sin(2t)/2+l/2) ^^ 

1 + ij27r/2 Jo 

I ^/^ f ''g-iiit-{l+ii27r/2)(sm{2t)/2-l/2)^^ 

1 + ii2^/2 Jo 



Again, we approximate these integrals with the fourth-order convergent composite Simpson's 
rule to get accurate approximations for the Fourier coefficients of q^. 
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Remark 5.1. Of course, a similar integration-by-parts trick with respect to xi would still work 
if 53 depends in a more complicated way on X2 ■ This shows that in principle the Fourier co- 
efficients of contrasts that vary smoothly in one variable can be computed by approximating 
one- dimensional integrals of smooth functions. 

Finally, we define the contrast function 54 plotted in Figure IHd) - a contrast function that 
varies smoothly in a 27r-periodic rectangle-shaped structure with support D, D = (—2.5, 2.5) x 
(-0.75,0.75). In detail, 

§4(2;) = 2cos(xi) (x2 + 0.75) for X = (3;i,X2) G -D 

and 54 (x) = for points outside of D. The Fourier coefficients of 54 can be explicitly computed 
using integration-by-parts techniques we already used above. Omitting technical details, the 
result is that 

«U) = 7S= fo^ •? = Ui'Ja) G Z , 

where 



Svr 



A[ji 



'sin(5ji)[(2cos(10)+l)/ii-8/i3]-4cos(5ii)sin(10)/i2 , ^ „ > .„ , „. 
^ YZi/Jf JlGZ\|0,±2}, 

sin(20)/4 + sin(lO) + 5 ji = ±2, 

[sin(10)/2 + 5 ji = 0. 



Bi,) = \^ exp(-37rijV4) - j^ sin(3^jV4) n + 0, 
^ ^^ l9/2 J2 = 0. 



Remark 5.2. The last example shows that Fourier coefficients of contrasts of the form 
q{x) = fi{xi)f2{x2) can be computed (semi-) analytically if fi^2 o'^e trigonometric functions, 
exponentials, or polynomials. The last example features a linear function f2{x2) = X2 + 0.75, 
however, higher- degree polynomials could be treated as well using additional integrations by parts 
reducing the polynomial degree. 

Since explicit analytic solutions for plane wave scattering problems involving the contrast 
functions 51,2,3,4 are not known, we check convergence rates for these structures by computing 
a reference solution for very large discretization parameter N. For all examples below, this 
reference solution is computed for N = 3072 using GMRES with a relative residual reduction 
factor of 10~®. The angle of the incident plane wave is always chosen as 9 = 7r/4 and the 
wave number always equals k = tt/2. We check the convergence rates from Proposition 14.11 
by computing scattered fields for discretization parameter N = 2", n = 4, . . . , 9. As above, 
the GMRES algorithm is stopped when the relative residual reduction factor is less than 10~^. 
Figure [5] shows that the convergence order of the method in the energy norm Hp^,,. is in good 
agreement with the statement of Proposition 14.11 Further, for all test cases, the rates of the 
error measured in -ffpcr and in L^ are in good agreement with the statement of Remark 14.21 
Computation times and the number of iterations of the GMRES algorithm corresponding to 
the numerical experiments illustrated in Figure [5] are shown in Table [2j 
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Figure 5: Test for the convergence rate of the trigonometric Galerkin discretization for the 
different structures presented in Figure |H The plots show the relative error in the H^ j.-norm 
between the approximate solution (N = 2", n = 4, . . . , 9) and the reference solution (N = 3072), 
plotted against the discretization parameter A^. Circles, kites, triangles correspond to s = 1, 
s = 0.5 and s = 0, respectively. The continuous line and the dotted lines indicate the expected 
convergence orders 0.5 and 1, respectively, (a) Results for the kite-shaped contrast qi from 
Figure m^a). (b) Results for the piecewise constant contrast (72 from Figure IHb). (c) Results 
for the contrast q^, that varies smoothly within a sinusoidal strip from Figure [D^c). (d) Results 
for the contrast q^ that varies smoothly within a rectangle from Figure [D^d). 



20 



N 


64 


128 


256 


512 


Time(s) for qi (Figure H^a)) 


1.4 


7 


44 


295 


Time(s) for g2 (Figure l^b)) 


1.7 


5 


16 


47 


Time(s) for qs (Figure Hl^c)) 


1.6 


7 


39 


184 


Time(s) for q^ (Figure IUd)) 


0.4 


2 


8 


39 


t] GMRES iterations for qi (Figure Hl^a)) 


10 


11 


11 


11 


tt GMRES iterations for gs (Figure ifb)) 


12 


12 


12 


12 


tt GMRES iterations for q^ (Figure ijc)) 


6 


6 


6 


6 


tt GMRES iterations for q^ (Figure H^c)) 


9 


10 


10 


10 



Table 2: Computation times and number of GMRES iterations for the computation of the error 
curves shown in Figure El Computing the reference solutions took roughly 1 hour for q2 and ^4 , 
10 hours for ^3 and 16 hours for qi. 

The last computational experiment illustrates the convergence of the trigonometric Galerkin 
technique using an error indicator resulting from energy conservation. Recall the Rayleigh 
coefficients u- of the scattered field from @. For the incident plane wave u* with incident 
angle 9, we define similar coefficients by u*- = J^^u^{xi, —h) exp{—iajXi) dxi for j G Z. Then 
Green's formula applied to equation ([TJ together with the Rayleigh expansion condition shows 
that 

/3j(l"/l' + l«7+4l')=/'o. (24) 






The sums 



i:fc2>/32 i:fc2>/32 

correspond to transmitted and reflected wave energies. In the following experiment, we compute 
the function 

9 ^ \1 - EtrM - Ercm\ (25) 

for many angles 9 to obtain an error indicator for the numerical accuracy of the integral equation 
solver in dependence on the angle of the incident field. This angle, 9, is sampled at 200 points 
uniformly distributed in the interval [0.2, 1.2]. The wave number k equals 2.5. To compute the 
energy curves shown in Figure E|^a) the scattered field is approximated in 7n where N = 2^ = 
256. The relative residual reduction factor for the GMRES iteration is in this experiment always 
chosen as 10~^. With this choice, the computation time for solving for one fixed incident angle 9 
is about 8 seconds. In Figure [6l^b) we check the error indicator of energy conservation from (|25p 
for different discretization parameters A^. This plot shows that the error of the computed 

Rayleigh coefficients corresponding to propagating modes converges with order 1, exactly as 

1/2 
the error of the solutions in Hpcr ■ This seems natural since, first, the Rayleigh coefficients are 

obtained from the numerical solution un by integration over the line Tp = {—it, tt) x {p} and, 

second, the trace theorem states that the mapping n^r 1— )• n^rlp is bounded from H^^j.{nfi) 
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into L'^iTp) for s > 1/2. The plot in Figure El^b) further shows a shght instabihty around a 
Wood's anomaly at the angle 6 = arccos(l — 1/(2.5)) w 0.927, as it is going to be expected from 
Remark 12.11 (The sampling points naturally avoid the exact value of the angle corresponding 
to this Wood's anomaly.) 
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Figure 6: (a) Reflected energy curves (low dashed line) and transmitted energy curves (contin- 
uous line) plotted against the angle 6 of the incident plane wave u*. The two curves sum up to 
one, as they should due to (^^. (b) The error criterion ([25|) plotted for different discretization 
parameters N = 2", n = 5, . . . , 10, versus the angle 9 of the incident plane wave. (The order of 
the curves from top to bottom corresponds to the increasing discretization parameter N.) 
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